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Abstract 

A  water  and  thermal  management  model  for  a  Ballard  PEM  fuel  cell  stack  was  developed  to  investigate  its  performance.  A  general  calculation 
methodology  was  developed  to  implement  this  model.  Knowing  a  set  of  gas  feeding  conditions  (i.e.,  pressure,  temperature,  flow  rate)  and  stack 
physical  conditions  (i.e.,  channel  geometry,  heat  transfer  coefficients,  operating  current),  the  model  could  provide  information  regarding  the 
reaction  products  (i.e.,  water  and  heat),  stack  power,  stack  temperature,  and  system  efficiency,  thereby  assisting  the  designer  in  achieving  the 
best  thermal  and  water  management.  Furthermore,  if  the  stack  undergoes  a  perturbation,  such  as  the  initial  start-up,  quick  change  in  current, 
or  a  shutdown,  the  model  could  predict  the  dynamic  information  regarding  stack  temperature,  cell  voltage,  and  power  as  a  function  of  time. 
©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

A  proton  exchange  membrane  (PEM)  fuel  cell  is  an  elec¬ 
trochemical  device  where  the  energy  of  a  chemical  reaction 
is  converted  directly  into  electricity,  by  combining  hydrogen 
fuel  with  oxygen  from  air  [1].  Water  and  heat  are  the  only  by¬ 
products  if  hydrogen  is  used  as  the  fuel  source  for  PEM  fuel 
cell.  Most  of  the  current  research  and  development  efforts  fo¬ 
cus  on  PEM  fuel  cells  due  to  their  capability  of  higher  power 
density  and  faster  start-up  than  other  fuel  cells  [2-6].  Usu¬ 
ally  PEM  fuel  cells  could  be  operated  at  a  temperature  lower 
than  100  °C,  thus  faster  start-up  and  immediate  response  to 
changes  in  the  demand  for  power  could  be  realized. 

Water  and  thermal  management  has  become  one  of  the 
key  technical  challenges  that  must  be  resolved  in  order  for 
the  PEM  fuel  cell  technology  to  be  feasible  for  transportation 
applications  [7,8],  although,  over  the  last  decade,  significant 
progress  has  been  made  in  the  field  of  PEM  fuel  cell  stack 
development  [9-1 1].  Proper  water  and  thermal  management 
is  essential  for  optimizing  the  performance  of  a  fuel  cell  stack. 
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In  automotive  applications,  there  are  many  different  road 
conditions  involved  and  therefore  the  knowledge  on  the  PEM 
fuel  cell  stack  in  terms  of  steady  and  transient  behavior  (e.g., 
acceleration,  deceleration)  becomes  very  important.  In  an  au¬ 
tomotive  fuel  cell  stack,  water  and  thermal  management  on 
this  steady  and  transient  behavior  is  associated  with  many  pa¬ 
rameters  that  affect  the  design  and  performance  of  PEM  fuel 
cell.  In  order  to  understand  the  relative  importance  of  the 
parameters  and  their  interaction,  an  investigation  on  these 
parameters  is  required  [12].  Mathematical  modeling,  a  con¬ 
venient  and  powerful  technique,  is  therefore  well  suited  for 
this  task.  The  numerical  modeling  could  be  employed  to  sig¬ 
nificantly  reduce  the  time  and  cost  associated  with  the  PEM 
fuel  cell  stack  development. 

To  date,  most  of  the  work  done  in  terms  of  PEM  fuel 
cell  modeling  has  focused  on  the  electrochemical  and  diffu¬ 
sion  processes  of  individual  fuel  cell  (also  called  unit  cell). 
Some  noteworthy  early  examples  include  Dunbar  and  Gag- 
gioli  [13],  Springer  et  al.  [14],  Verbrugge  and  Hill  [15], 
Bernardi  and  Verbrugge  [16,17],  Fuller  and  Newman  [18], 
Nguyen  and  White  [19]  and  Kim  et  al.  [20].  University  of 
Victoria  and  University  of  Waterloo  [21-25]  have  been  con¬ 
ducting  the  fuel  cell  modeling  for  many  years  and  have  made 
very  impressive  progress  on  the  unit  cell  modeling. 
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Nomenclature 

a  water  activity 

A  area  (m2) 

c  water  concentration  in  the  membrane 

(molm-3) 

Cp  average  heat  capacity  (J  kg- 1  K- 1 ) 

d  channel  height  (m) 

Dm  water  diffusion  coefficient  (m2  s-1) 

E  thermodynamic  potential  (V) 

/  fraction  coefficient  of  channel 

F  Faraday’s  constant  (F  =  96,485) 

h  convective  heat  transfer  coefficient 

(W  m  2  K-1) 

H  change  in  enthalpy  (J  mol-1) 

AH  heat  of  reaction  (J  mol- 1 ) 

i  operating  current  (A) 

I  current  density  (A  cm-2) 

/o  exchange  current  density  (A  cm-2) 

kp  water  hydraulic  permeability  in  membrane 
(m2  s-1) 

K  thermal  conductivity  (W  s-1  m-1  K-1) 

F  channel  length  (m) 

M  mass  of  the  fuel  cell  (kg) 

n  number  of  cells  in  the  stack 

rid  electro-osmotic  drag  coefficient 

nQ  mole  number  of  electrons  per  unit  current  per 
unit  time 

N  molar  flow  rate  (mol  s-1);  channel  number 
P  power  output  (W) 

PEM  proton  exchange  membrane 

q  energy  (W) 

R  universal  gas  constant  (=8.314  J  mol-1  K-1) 
RH  relative  humidity 

t  thickness  (m) 

T  temperature  (K) 

V  output  voltage  (V);  velocity  (ms-1) 

Xi  mole  fraction  of  species  i 

Greek  letters 

a  excess  coefficient 

v)  over-voltage  (V) 

X  water  content  of  membrane 

/j  water  viscosity  (Pa  s) 

0  relative  water  content 

Subscripts 

0  standard  condition 

a  anode 

act  activation 

c  cathode 

cell  proton  exchange  membrane  fuel  cell 
cons  consumed 

conv  convection  flux 


diff 

diffusion  flux 

elec 

electrical 

g 

gas 

hum 

humidification 

in 

in 

inlet 

flow  inlet  stack  channel 

int 

internal 

1 

liquid 

loss 

loss 

m 

membrane 

mass 

mass  transfer  and/or  mass  consumption 

ohmic 

ohmic 

out 

out 

outlet 

flow  outlet  stack  channel 

prod 

product 

room 

ambient  condition 

rxn 

reaction 

sens 

sensible 

stack 

fuel  cell  stack 

theo 

theoretical 

trans 

water  transfer  across  membrane 

w 

water 

Superscripts 

a vg 

average  value 

channel 

stack  flow  channel 

dry 

dry  gas  condition 

new 

current  value  in  iterative  calculation 

old 

previous  value  in  iterative  calculation 

sat 

saturation  condition 

* 

at  the  catalyst  interface 

The  models  mentioned  above  mainly  emphasized  on  un¬ 
derstanding  and  improving  the  kinetic  processes  that  oc¬ 
curred  in  fuel  cell,  aiming  at  improving  individual  fuel  cell 
performance.  The  researchers  built  their  models  based  on 
electrochemical  theories,  electrode  kinetics  and  experimen¬ 
tal  data. 

As  mentioned  by  Costamagna  and  Srinivasa  [26],  until  the 
year  2000,  no  detailed  results  of  the  modeling  analyses  of  the 
performance  characteristics  of  the  electrochemical  cell  stack 
and  the  PEMFC  power  plant  had  appeared  in  the  literature. 
Models  of  fuel  cell  stacks  have  been  and  are  being  conducted 
by  some  fuel  cell  companies  and  such  development  remains 
in  proprietary. 

Texas  A&M  University  [27,28]  made  very  good  contri¬ 
bution  to  the  fuel  cell  stack  modeling.  However,  their  model 
only  focused  on  fuel  cell  stack  and  the  model  did  not  con¬ 
sider  two-phase  flow  and  liquid  water  was  not  considered.  In 
real  fuel  cell  processes,  both  liquid  water  and  vapor  are  very 
important  factors  that  have  to  be  resolved  properly  in  order 
to  have  stable  fuel  cell  operation. 
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Some  thermal  models  of  PEM  fuel  cell  stacks  could  be 
found  in  the  literature  [29-32].  These  models  typically  treat 
the  stack  as  a  process  unit  and  develop  models  based  on  elec¬ 
trochemical  performance,  and  the  physical  characteristics  of 
the  inlet  and  outlet  flows.  The  computations  of  these  models 
are  usually  too  involved  to  be  employed  in  a  comprehensive 
model  of  a  PEM  fuel  cell  stack.  A  need  exists  for  a  technique 
that  could  be  used  to  determine  PEM  fuel  cell  stack  thermal 
performance  without  requiring  significant  amount  of  com¬ 
putations.  Some  excellent  studies  on  these  topics  have  been 
conducted  by  a  group  of  scientists  in  Royal  Military  College 
of  Canada  [33-36].  In  [37]  by  Yu  and  Zhou,  an  improved 
model  was  built  to  consider  the  inlet  water  vapor  effects. 

To  the  author’s  knowledge,  the  models  mentioned  above 
have  not  included  the  liquid  water  effects,  especially  the  inlet 
water  (liquid  and  vapor)  effects  that  could  play  a  very  im¬ 
portant  role  in  the  PEM  fuel  cell  performance.  Therefore,  in 
the  present  study,  a  two-phase  model  with  phase  change  was 
built  to  meet  this  challenge. 


2.  Model  development 


2.7.  Steady-state  electrochemical  model 


The  steady-state  electrochemical  model  could  be  used  to 
predict  stack  voltage  output.  The  cell  voltage  was  defined  in 
terms  of  the  following  three  terms  [33]:  the  thermodynamic 
potential  E ,  the  activation  over-voltage  77act,  and  the  ohmic 
over-voltage  ^hmic : 

kcell  =  E  —  T]act  —  7?ohmic  (2) 

where 

E  =  1.229  -  0.85  x  10“3  x  <Tstack  -  298.15) 

+  4.3085  x  1(T5  x  Tstack  x  Pn(P&2)  +  0.5  x  ln(^2)] 

(3) 


Here  the  friction  effect  (pressure  drop)  was  introduced  to 
modify  the  original  model  [33]  by  averaging  the  inlet  and 
outlet  partial  pressures: 


PU2  2^H2,in  +  PH2,out) 
Po2  =  fp*02,in  +  Po2,  out) 


(4) 

(5) 


For  modeling  purposes,  the  following  assumptions  were 
made  in  the  present  study: 


The  activation  overpotential  and  ohmic  overpotential  could 
be  calculated  as  follows  [19]: 


(1)  The  product  water  generated  at  the  cathode  is  assumed 
to  be  in  liquid  state. 

(2)  The  water  condensation/evaporation  rate  is  not  consid¬ 
ered.  Instead,  the  equilibrium  between  the  water  vapor 
and  liquid  is  always  assumed. 

(3)  Ideal  gas  law  was  employed  for  gaseous  species. 

(4)  Stack  temperature  is  uniform  due  to  high  thermal  con¬ 
ductivity. 

(5)  Water  transport  in  and  out  of  the  electrodes  was  in  the 
form  of  vapor. 

(6)  The  electrode  layers  were  “ultra- thin”,  so  that  gas  trans¬ 
port  resistance  through  the  electrode  porous  layer  could 
be  neglected. 

(7)  The  liquid  water  was  assumed  to  exist  at  the  surface  of 
the  channels,  and  the  volume  to  be  negligible. 

(8)  For  the  pressure  drop  calculation,  the  liquid  water  was 
neglected.  The  entrance  and  exit  losses  were  neglected, 
which  were  too  small  compared  with  the  overall  pressure 
drop. 

In  order  to  describe  both  cases  either  with  or  without  phase 
change,  a  parameter  0,  relative  water  content,  was  defined  as 
follows: 

Total  mole  number  of  water  (vapor  +  liquid) 
Maximum  possible  mole  number  of  water  vapor 

According  to  assumption  (2),  when  0  <  1,  it  is  exactly  the 
same  as  relative  humidity  and  there  is  no  liquid  water;  while 
0  >  1  means  there  is  liquid  water  and  0  is  no  longer  equivalent 
to  the  relative  humidity. 


7?  act  — 


7?(273.15  +  rstack) 
0.5  F 


In 


0  ohmic  — 


Itm 


(6) 

(7) 


where  rstack  is  the  stack  temperature  (K),  i  the  operating  cur- 
rent  (A),  p  the  partial  pressure  on  the  catalyst  interfaces  cor¬ 
responding  to  concentration  of  feeding  gas,  and  /  the  current 
density. 


2.2.  Steady-state  thermal  model 

A  steady- state  thermal  model  was  established  based  on 
the  balance  of  mass  and  energy  about  fuel  cell  stack.  Fig.  1 
shows  a  schematic  of  the  inlet  and  outlet  streams  in  a  typical 
PEM  fuel  cell  system.  Hydrogen,  air,  and  cooling  water  are 
independent  streams.  Energy  balance  about  the  fuel  cell  stack 
was  performed  to  calculate  various  energy  terms  associated 
with  fuel  cell  operation: 


^theo  —  ^elec  T  ^sens  T"  ^latent  T  ^loss  (8) 

where  gtheo  is  the  theoretical  energy  produced  by  the  fuel  cell 
reaction,  qsens  the  sensible  heat  calculated  for  each  of  the  fuel 
cell  streams  (anode,  cathode,  and  water  coolant),  giatent  the 
total  latent  heat  of  the  water  vaporization  (condensation)  for 
anode  and  cathode  streams,  geiec  the  electrical  energy  output, 
and  gioss  the  heat  loss  from  the  surface  of  the  stack.  Compar¬ 
ing  (8)  with  the  model  used  in  [33] ,  the  model  developed  in  the 
present  work  included  the  two-phase  effect  (phase  change). 
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fltheo 

V 


N  H2,a,in 
N  C02„a,in 
N  w,g,a,in 

N  w,l,a,in 

T  P 

1  a, m,  r  a,in 


N  w,in  T  w  jn 


N  02,c,in 
N  N2,c,in 
N  w,g,c,in 
N  w,l,c,in 


1  1  i  1 

flsens  Q  latent  flloss  Qelec 


FUEL  CELL 
STACK 

Tstack 


N  H2,a,out 
N  C02,a,out 
N  w,g,a,out 
N  w,l,a,out 
T  a, out ,  P  a, out 


N  w„out  T  w,out 


N  02,  c,  out 
N  N2,c,out 
N  w ,g,c,out 
N  w,l,c,out 
Tc,OUt?  Pc, out 


Fig.  1.  Schematic  of  the  inlet  and  outlet  streams  in  a  fuel  cell  system  iden¬ 
tifying  various  parameters  including  flow  rates  and  energy  terms. 


2.2.1.  Energy  equations 

Theoretical  energy  from  the  electrochemical  reaction  in 
PEM  fuel  cell  was  calculated  through  the  product  of  reaction 
energy  and  molar  flow  rate  of  consumed  hydrogen: 

^theo  =  ^F[2,cons  A  Hrxn  (9) 

The  electrical  power  generated  by  the  PEM  fuel  cell  stack 
with  n  single  cells  was  evaluated  as 

^elec  —  ft  Ecelfl  (10) 

The  sensible  heat  through  anode  stream  was  considered  for  all 
the  possible  species  in  anode  (i.e.,  water  vapor,  liquid  water, 
carbon  dioxide  from  reformate)  as  follows: 


^sens,a  =  Wh2,  a,out^p,H2,g(^a,out  '  -T0) 

+  ^W,g,  a,  out  Cp,  H2O,  g  (  T&,  out  To) 
+  Vco2,  a,  out  C  CO2 ,  g  (  ^a,  out  To) 
“I-  ^W,l,a,0ut^/?,W,l,out(Ta,out  —  ^b) 

—  ^w,l,a,in^/?,w,l,in(^a,in  —  ^b) 

—  ^w,g,a,inC/?,H20,g(^a,in  —  7q) 


por  and  water  liquid  were  assumed  to  be  in  equilibrium  all  the 
time,  i.e.,  the  condensation/evaporation  process  was  assumed 
to  be  so  fast  that  there  is  no  finite  condensation/evaporation 
rate;  also  the  water  transfer  across  the  membrane  was  as¬ 
sumed  in  vapor  form,  see  details  for  this  assumption  in  [19]): 


(12) 


The  sensible  heat  in  cathode  was  considered  in  a  similar  way 
to  that  in  anode  except  the  species  are  different  from  those 
in  anode.  In  cathode  the  species  include  oxygen,  nitrogen, 
water  vapor,  water  liquid,  as  shown  in  (13): 


^sens,c  —  No2,c  ,out  C  p,  02,g  (Tc  ,out  -To) 

+  ^W,g,  c,out  CPi  FUO,g(^c,out  To) 

+  An2,  c,out^p,N2,g(^c,out  To) 

T  ^w,l,c,out^p,w,l,out(^c,out  —  To) 

—  ^w,l,c,in^/?,w,l,in(^c,in  —  Tq) 
^w,g,c,inCp,H20,g(^c,in  —  Tq) 

—  ^02,c,in^p,02,g(^c,in  —  ^b) 

—  ^N2,c,inCp,N2,g(^c,in  —  ^b) 


(13) 


The  latent  heat  in  cathode  is  somehow  complicated  due  to 
the  water  generation,  water  phase  change  and  transfer  across 
the  membrane.  The  basic  rule  here  is  to  figure  out  the  molar 
flow  rate  of  the  water  vapor  that  is  involved  in  phase  change. 
Details  are  as  follows. 

For  latent  heat  in  cathode,  if  Nw,l,c,in  >  (Nw,g,c,out  - 
A^rans  —  ^w,g,c,in)»  i.e.,  the  amount  of  liquid  water  carried 
from  the  cathode  inlet  is  big  enough  for  phase  change,  then 
we  have 


^latent, c  —  (^w,g,c,out  ^trans  ^w,g,c,in)^vaporizationl,cl 


(14) 


Otherwise,  the  liquid  water  carried  from  the  inlet  must  be 
evaporated  and  some  of  product  water  must  be  evaporated 
too,  so  we  have 


^latent,  c  —  ^w,l,c,  in  ^vaporization,  cl  T  (^w,g,c,out  ^trans 

—  ^w,g,c,in  —  N w,  1,  c,  in)  ^vaporization,  c2  (13) 


where 

^vaporization  =  45070  -  41.97  +  3.44  X  10_3r2 

+  2.54  X  10“67’3  -  8.98  x  10“107’4  (16) 


—  ^H2,a,inCp,H2,g(Ain  —  To) 

Vco2,a,inCpiC02,g(T,in  T))  (U) 

The  latent  heat  through  the  anode  was  included  through  track¬ 
ing  down  the  phase  change  (in  the  present  paper,  the  water  va- 


and  subscripts  cl  and  c2  represent  the  different  state  (thus 
different  temperature)  for  the  evaporation  of  water  that  are 
from  different  origin,  e.g.,  either  from  inlet  stream  or  electro¬ 
chemical  product.The  sensible  heat  in  water  coolant  stream 
was  calculated  by  use  of  the  following  formula: 
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^sens,w  —  N w,inC/?,w,l(^w,out  A)) 

^w,inC/7,w,l(^w,m  —  Tq) 


(17) 


Then  the  sensible  and  latent  heat  were  summed  and  the  heat 
loss  from  the  stack  to  the  ambient  was  calculated  based  on 
(8): 


^sens  =  ^sens,a  T  ^sens,c  T  ^sens,w 

^latent  =  ^latent,  a  T  ^latent,  c 

^loss  =  ^theo  —  ^elec  —  ^sens  —  ^latent 


(18) 

(19) 

(20) 


2.2.2.  Flow  rates 

The  saturation  pressure  (atm)  was  calculated  based  on  the 
following  equation  [19]: 

sat  i  q— 2.1794+0.029537— 9. 1837x  10-5r2  +  1.4454x  10-7r3 
Pw,g 

(21) 


The  molar  flow  rate  for  hydrogen  in  anode  and  air  in  cathode 
on  dry  condition  at  each  inlet  can  be  evaluated  according 
to  the  operating  current  and  excess  coefficient  [37]  on  each 
stream  inlet: 


1 

^a,H2,in,dry,0  —  777  Ioin2ne 

2pH2 

1 

Nq,  air,  in,  dry,  0  —  ~T7 

4P02 


(22) 

(23) 


where  nQ  =  1.0365  x  10  5  mol/A  s  is  the  molar  flow  rate  of 
electrons  for  generating  1  A  electricity,  a  the  excess  coeffi¬ 
cient,  i.e.,  the  ratio  of  the  actual  amount  supplied  to  the  the¬ 
oretical  amount  needed,  and  /3  the  molar  fraction  of  oxygen 
in  air  stream  at  cathode  inlet. 

The  equations  of  flow  rates  were  proposed  to  account  for 
the  inlet  water  (liquid  +  vapor),  as  listed  below. 

The  maximum  water  vapor  carried  from  the  anode  inlet 
was  evaluated  as 


A 


w,g,a,m,max 


—  (^H2,a,in  T  AcQ2,a,in) 


P 


sat 

w,g,a,in 


P a, in  -  P 


sat 

w,g,a,in 


(24) 


Then  the  amount  of  water  vapor  and  water  liquid  at  the  anode 
inlet  were  calculated  as  below. 

If  Aw,a,in  2  Aw,g,a,in,max?  then  we  have 


P^w,g,a,in  —  AW5g?a5in5max, 
^w,l,a,in  =  ^w,a,in  —  ^w,g,a,in 


(25) 


If  Aw?a5in  <  AW  g?a,in,max  ?  then  we  have 


^w,g,a,in  —  (^H2,a,in  T  ^C02,a,in) 


Pa, in  Pw,g,a,in0a3n 


^w,l,a,in  —  0 


(26) 


The  maximum  amount  of  water  vapor  at  anode  outlet  was 
calculated  as  follows: 


A 


w,g,a,out,max 


—  (AH2,a,out  T  AcQ2,a,out) 


P 


sat 

w,g,a,out 


sat 


(27) 


Pa, out  Pw,g,a,out 


Then  the  amount  of  water  vapor  and  water  liquid  in  the  anode 
outlet  were  evaluated  as  below. 

If  Aw?a,in  —  Atrans  ^  ^w,g,a,out,max?  then  we  have 


^w,g,a,out  —  ^w,g,a,out,max? 

^w,l,a,out  =  ^w,a,in  —  Atrans  —  ^w,g,a,out 

If  AW5a?in  —  A^ans  <  ^w,g,a,out,max?  then  we  have 


(28) 


^w,g,a,out  —  ^w,a,in 


trans 


A 


w,l,a,out 


=  0 


(29) 


For  cathode  inlet,  the  maximum  water  vapor  carried  from  the 
cathode  inlet  was  evaluated  as 


A 


w,g,c,m,max 


—  (^C>2,c,in  T  AN2,c,in) 


P 


sat 

w,g,c,in 


Pc, in  Pw,g,c,in 


(30) 


Then  the  amount  of  water  vapor  and  water  liquid  in  the  cath¬ 
ode  inlet  were  evaluated  as  below. 

If  Aw  c  in  >  Aw?g  C  in?rnax?  then  we  have 


^w,g,c,in  —  ^w,g,c,in,max!> 
^w,l,c,in  =  ^w,c,in  —  ^w,g,c,in 


(31) 


If  Aw?c?in  <  AW  g?c?in?rnax?  then  we  have 


^w,g,c,in  —  (^02,c,in  T  ^N2,c,in) 


P 


sat 
w,g,c,m 


in  0c,  in 


Pc, in  Pw,g,c,in0cdn 


A 


w,l,c,in 


=  0 


(32) 


In  cathode  stream,  the  water  was  produced  and  the  product 
water  was  assumed  to  be  liquid  in  the  present  study.  It  was 
evaluated  as 


^w,l,prod  —  An2,cons  —  ^H2,a,in  ^H2,a,out  (33) 

For  cathode  outlet,  the  maximum  water  vapor  carried  from 
the  cathode  outlet  was  evaluated  as 

_sat 

^w,g,c,out,max  —  (^C>2, c, out  T  ^N2,c,out)  ^ 

Pc, out  —  P  w,g,c,out 

(34) 

Then  the  amount  of  water  vapor  and  water  liquid  in  the  cath¬ 
ode  outlet  were  evaluated  as  below. 

If  (Aw?C5in  +  Awj?prod  +  Atrans)  2  AW5g?C50Ut?max>  then  we 
have 


^w,g,c,out  —  ^w,g,c,out,max> 

^w,l,c,out  =  ^w,c,in  T  ^w,l,prod  T  ^trans  —  ^w,g,c,out  (35) 

If  (Aw?C5in  +  Awj?prod  +  Atrans)  <  ^w,g,c,out,max?  then  we 
have 
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^w,g,c,out  —  ^w,c,in  T  ^w,l,prod  T  A^rans? 
^w,l,c,out  =  0 


(36) 


The  average  heat  transfer  coefficient  for  the  stack  may  be 
estimated  using  the  average  heat  loss  from  the  surface  of  the 
fuel  cell  stack.  Similarly,  the  increase  in  sensible  and  latent 
heat  terms  could  also  be  linked  to  heat  transfer  coefficients,  hj , 
from  the  stack  to  the  fluid  j,  where  j  =  anode,  cathode,  or  water 
stream.  Once  heat  transfer  coefficients  h ,  heat  exchange  area 
and  sensible,  latent  heat  terms  are  known,  the  temperature 
of  stack  and  outlet  flows  could  be  estimated  by  using  the 
following  equations: 


Tstack  — 


^loss 


(/z  A)  stack 


+  rf 


room 


(37) 


Z 


a,  out 


Tstack 


4s&n,a.  T  ^latent, a  T  ^mass,a 

(hA\ 


(38) 


all  the  terms  on  the  right  side  of  the  Eq.  (44),  it  could  be  used 
as  a  basis  of  a  finite-difference  calculation  using 


'rmew 
1  stack 


yOld 
1  stack 


+ 


Attack 
d  t 


At 


(45) 


2.4.  Pressure  drop 


Pressure  drop  along  the  channels  could  be  calculated  by 
using  average  gas  velocity,  which  is  the  mean  value  of  inlet 
and  outlet  velocity  of  each  stream.  Ignoring  the  volume  of 
liquid  water,  the  local  velocity  V  (m  s-1)  was  determined  by 
gas  molar  flow  rate  (mols-1),  local  pressure,  temperature, 
cross-section  area  of  channel  Ac,  and  number  of  channels 
(Nch): 


N  x  22.4  x  1CT3  po  T 
^4-cA^ch  P  Zq 


(46) 


where  the  gas  molar  flow  rate  could  be  determined  for  each 
stream  as  follows. 


Z 


c,out 


=  2 


Tstack 


T  ^latent, a  ^mass^ 

(hA)c 


-Z 


c,m 


(39) 


7w,OUt 


Tstack 


^sei^w 

(hA)w_ 


-  Z 


w,m 


(40) 


where  the  energy  change  due  to  mass  transfer  and  mass  con¬ 
sumption  (including  the  sensible  energy  carried  by  the  water 
transfer  across  the  membrane,  the  sensible  energy  carried  by 
hydrogen/oxygen  consumed)  was  evaluated  as  follows: 


^mass^ 


A^trans6^p,H20,g(7stack  Zq) 

+  Nu2  ,con  Cp,  H2,g(7stack  7q) 


(41) 


at  anode  inlet: 


A^  —  (A^H2,a,in  T  A^C02,a,in)  I  1  T 


^sat  p  T_r 

TV, g,a,in  Rna,m 
Pa, in  -  Pw,g,a,in  RHa,i 


,in 


(47) 


at  anode  outlet: 


N  —  (A^H2,a,out  T  A^CQ2,a,out) 


.sat 


X  1  + 


Pw,g,a,out  RHa,out 


Pa, out  Pw,g,a,out  RHaoUt 


.sat 


(48) 


^mass^  —  A^trans^7p,H20,g(7stack  Zq) 

+  ^Vh2  ,con  Op,  H2  0,1  (Tstack  Zq) 

—  N o2,conC p ,02( Tstack  ~  Zq)  (42) 

2.3.  Transient  model 


In  transient  state,  an  additional  accumulation  term  should 
be  considered,  therefore: 


m  stack  Cp,  stack 


fiT'^ck 
d  t 


—  ^theo  ^elec  ^sens  ^latent  4\oss 


(43) 


where  m  is  the  total  mass  of  the  fuel  cell  stack,  C  the  aver¬ 
age  specific  heat  of  the  stack,  and  drstackAfi  the  temperature 
change  with  respect  to  time.  From  Eq.  (43),  we  have 


d  Ts  tack 


^theo  ^elec  ^sens  ^latent  ^loss 


d  t 


m  stack  Cp,  stack 


(44) 


In  the  calculations  presented,  an  average  value  of  35  kJ  K  1 
was  used  for  m  stack1 Cp,  stack  of  Ballard  Mark  V  stack.  Knowing 


at  cathode  inlet: 


N  —  (A^Q2,c,in  T  A^N2,c,in) 


P 


sat 

w,g,c,in 


RH 


c,m 


Pc, in  -  P 


sat 

w,g,c,in 


RH 


c,m 


at  cathode  outlet: 


N  —  (A^02,C,OUt  T  A^N2,C,OUt) 


X 


nsat 

Pw,g,c,out 


RH 


c,out 


Pc,  out  P 


sat 

w,g,c,out 


RH 


c,out 


(50) 


When  RH=  1,  the  largest  molar  flow  rate  for  each  stream 
is  obtained. Once  temperature  and  flow  rate  are  known,  the 
pressure  drop  along  the  channels  could  be  obtained  by  using 
(Darcy- Weisbach  equation  [38]): 


Ea  Pa  E 


a,m 


2 


2 


(51) 

(52) 


where  D  is  the  hydraulic  diameter. 
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2.5.  Water  transfer  across  membrane 

Water  transfer  across  membrane  is  the  sum  of  following 
three  terms  [19,39]: 

1 .  Electro-osmotic  drag  flux,  which  is  caused  by  hydrogen 
ion  drag. 

2.  Diffusion  flux,  which  is  caused  by  water  concentration 
gradient  between  anode  and  cathode. 

3.  Convection  flux,  which  caused  by  pressure  gradient. 

Therefore, 


Convection  flux  was  calculated  as  follows: 


_  kp  dpv 
iVconv,y  —  c  , 

M  dy 


IJ  d  y 


(61) 


where  kv,  /a,  d pv  and  c  are  the  hydraulic  permeability  of  wa¬ 
ter  in  membrane,  water  viscosity,  partial  pressure  difference 
between  anode  and  cathode,  and  concentration  of  water  in 
membrane. 


3.  Solution  methodology 


N trans  —  2/ drag  T  51 d i  jf  +  N COnv 


(53)  3.1.  Steady -state  models 


The  electro-osmotic  drag  flux  could  be  calculated  by  [19,39] : 


51  drag  —  nd 


m 

F 


(54) 


n  d 


(55) 


X  =  0.043  +  17.81a  —  39.85a2  +  36.0a3  at  (a  <  1) 

=  14.0  +  1.4(a  —  \)at  (3  >  a  >  1)  =  16.8  at  (a  >  3) 

(56) 


Step  1 .  Start  with  a  guess  or  estimate  for  the  values  of  Tstack, 
7w,out?  Ta,out?  ^c,out  and  pout. 

Step  2.  From  these  guessed  values,  calculate  tentative  values 
of  related  energy  terms  and  Vceii. 

Step  3.  Use  those  tentative  energy  values  and  heat  transfer 
coefficients  to  get  new  values  of  T  and  p. 

Step  4.  With  these  p ,  T  as  better  guesses,  return  to  step  2, 
repeat  the  process  until  further  repetitions  cease  producing 
any  significant  changes  in  these  values. 

Step  5.  These  final  values  of  T,  p  will  satisfy  energy  and  mass 
balance,  and  will  be  the  steady-state  result  of  the  stack. 
Step  6.  Other  related  values  of  parameter  can  be  calculated 
from  them. 


Pv  apor 
a  =  — — 

Psat 


(57)  32.  Unsteady -state  models 


where  n&  is  called  electro-osmotic  drag  coefficient,  I  the 
current  density,  a  the  water  vapor  activity  (ratio  of  the  wa¬ 
ter  vapor  pressure  and  the  saturation  pressure),  A  the  water 
content  of  membrane  that  is  related  with  water  vapor  activ- 
ity.Diffusion  drag  flux  is  decided  by  diffusion  coefficient  Dm, 
water  concentration  c  and  the  membrane  charge  concentra¬ 
tion  Cf  which  is  fixed  for  one  type  of  membrane  [19,39]. 


10  10  exp 

+  0.0264A 
10-10  exp 


2416  (sa)  -  (?) 

2  -0.00067U3)at(A 

2416  (555)  -  (?) 


(2.563  -0.33A 

>4) 

(-1.25A 


+  6.65)  at  (4  >  A  >  3) 


10  10  exp 


2416 


(2.05A 


The  time  step  At=  1  s  was  used  in  the  dynamic  calcula¬ 
tions,  thus  changes  of  all  the  parameters  could  be  traced  at 
each  second. 

Step  1.  Calculate  energy  term  and  Vceii  by  initial  input  val¬ 
ues.  Use  unsteady- state  thermal  model  equation  to  get  the 
value  of  drstack/dt  at  the  beginning  of  the  first  time  step. 
Step  2.  Calculate  rstack  value  at  the  end  of  the  first  time  step, 
guess  the  value  of  Twout,  Ta?ou t,  Tc  ou t  and p out¬ 
step  3.  Keep  fixed  value  of  rstack,  follow  steady-state  cal¬ 
culation  steps  to  find  all  parameter  values  at  the  end  of  the 
first  time  step. 

Step  4.  For  next  time  step,  go  to  step  1,  use  those  value  got 
from  step  3  as  the  initial  values,  and  repeat  the  process  until 
reach  the  end  of  time  period. 


4.  Results  and  discussions 


—  3.25)  at (3  >  A  >  2) 


(58) 


4.1.  Steady  cases 


c  =  Acf 


51 diff, 


dc 

d y 


dA 

d y 


In  the  calculations  presented  here,  unless  specified  other- 

(59)  wise,  RH  =  1  for  both  anode  and  cathode  inlet  stream.  Table  1 
shows  the  input  data  for  the  calculated  case  that  was  similar 

(60)  to  the  case  reported  by  Amphlett  et  al.  [36]  and  a  comparison 
has  been  discussed  by  Yu  and  Zhou  [37]. 
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Table  1 


Input  data  for  calculated  case  (/=  20  A) 


Parameter 

Value 

Aht2,a,in  (mol  s_1) 

0.0078 

Fa, in  (°C) 

23.5 

A, in  (Psig) 

35 

-V 02,c,in  (mol  s-1) 

0.004 

To, in  (°C) 

23.5 

Pc, in  (Psig) 

35 

AT, in  (mols-1) 

1.84 

Tv, in  (°C) 

23.5 

Toom  (°C) 

23.5 

Tell 

35  cells 

Fig.  2  shows  temperature  of  the  exits  at  anode,  cathode, 
and  water  coolant  with  respect  to  the  steady  operating  cur¬ 
rents  from  1  to  60  A.  It  could  be  seen  that  all  the  temper¬ 
atures  increased  with  the  increase  of  steady  operating  cur¬ 
rent  and  the  cathode  exit  temperature  was  higher  than  the 
stack  temperature,  anode  temperature  and  coolant  tempera¬ 
ture. 

Fig.  3  is  the  output  stack  power  at  different  steady  op¬ 
erating  currents  from  1  to  60  A.  The  power  output  almost 
increased  linearly  with  the  steady  operating  current. 

4.2.  Steady  cases  with  different  values  off  at  inlets 

In  Fig.  4,  the  anode  exit  temperature  at  0a,in  =  1.0  and  1.5 
was  plotted.  Here  0a,in>  relative  water  content  at  anode  inlet, 
represents  the  molar  ratio  between  total  amount  of  supplied 
water  (liquid  +  vapor)  at  anode  inlet  and  the  saturated  water 
vapor  carried  by  the  anode  inlet  stream.  When  0a?in  <  1 ,  the 
anode  outlet  temperature  did  not  vary  significantly  with  0a  in 
(Fig.  4a).  When  </>a,in  >  1 ,  liquid  water  would  mix  with  anode 
inlet  stream  and  thus  different  0a?in  values  would  have  an  ob¬ 
vious  effect  on  anode  outlet  temperature,  attributable  to  liquid 
water  vaporization  leading  to  anode  exit  temperature  reduc- 


Fig.  2.  Stack  temperature  and  exit  temperatures  at  anode,  cathode,  and  water 
coolant  with  steady  operating  currents  from  1  to  60  A. 


Fig.  3.  Stack  output  power  with  steady  operating  current  from  1  to  60  A. 


Fig.  4.  (a)  Anode  exit  temperature  at  0a,in  =  O.5  and  1.0.  (b)  Anode  exit 
temperature  at  different  0a,in  value  1.0  and  1.5. 
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Current  (A) 

Fig.  5.  Stack  temperature  and  voltage  with  different  steady  operating  current 
from  1  to  60  A. 

tion  (Fig.  4b).  Basically,  when  inlet  0a?in  value  increased,  the 
anode  exit  temperature  decreased. 

Fig.  5  shows  the  stack  temperature  and  voltage  with  differ¬ 
ent  operating  currents.  Output  voltage  decreased  when  cur¬ 
rent  was  increased,  attributable  to  a  higher  current  creating 
a  larger  ohmic  over-voltage  loss.  Fig.  6  shows  stack  volt¬ 
age  output  at  different  </>a,in  for  steady  operating  conditions. 
Electro-osmotic  drag  would  be  the  dominant  factor  affecting 
the  amount  of  water  transferred  across  the  membrane.  Water 
was  dragged  from  the  anode  to  cathode  side  resulting  in  dry 
gas  at  the  anode  side  which  would  reduce  membrane  conduc¬ 
tivity  and  subsequently  lower  the  stack  voltage.  Therefore,  in 
order  to  achieve  a  higher  voltage  output,  extra  humidification 
has  to  be  provided  to  the  gas  at  the  anode  side;  It  also  could 
be  seen  in  Fig.  6  that  voltage  output  increased  when  the  inlet 
0a, in  was  raised. 


Fig.  6.  Stack  voltage  with  different  0a,in  at  steady  operating  current  from  1 
to  60  A. 


Fig.  7.  Transient  exit  temperature  plots  of  the  start-up  process  for  the  oper¬ 
ating  current  at  30  A. 

4.3.  Unsteady  cases 

Figs.  7  and  8  show  the  start-up  characteristics  of  this  stack 
at  the  operating  current  of  30  A;  Fig.  7  shows  the  transient  exit 
temperature  plots  while  Fig.  8  shows  the  stack  temperature 
and  voltage. 

From  Figs.  7  and  8,  it  could  be  seen  that  the  stack  required 
about  30-40  min  to  reach  steady  state  with  the  operating  cur¬ 
rent  of  30  A.  In  the  first  20  min,  the  rate-of-exit  temperature 
increase  was  high  then  slowly  reduced  until  about  40  min 
when  the  stack  almost  reached  its  steady  operating  state. 

The  transient  response  of  the  stack  for  the  load- set-up  from 
30  A  (for  60  min)  to  50  A  (for  another  60  min)  are  shown  in 
Figs.  9  and  10.  In  general,  the  stack  required  about  40  min  to 
reach  its  steady  operating  state  after  the  load  was  changed. 
It  is  noteworthy  that  when  the  load  was  changed  from  30  to 
50  A,  the  immediate  exit  temperatures  of  cathode  and  water 


Fig.  8.  Transient  plot  for  the  stack  temperature  and  voltage  of  the  start-up 
process  for  the  operating  current  at  30  A. 
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Fig.  9.  Transient  plot  of  exit  temperatures  at  anode,  cathode  and  water 
coolant  during  the  load  set-up  from  30  to  50  A. 

coolant  decreased  because  the  amount  of  air  and  water  with 
lower  temperature  increased  from  the  inlets.  Furthermore, 
from  Fig.  10,  the  stack  voltage  at  operating  current  of  50  A 
was  lower  than  that  at  30  A,  due  to  the  increase  of  current 
that  would  create  larger  over- voltages  and  thus  smaller  cell 
and  stack  voltage. 

Fig.  11  shows  the  temperature  change  in  all  streams  as 
function  of  time  and  current  change  while  Fig.  12  gives  the 
stack  voltage  output  and  temperature  change  as  function  of 
time  and  load  change.  The  load  changed  in  each  20  min,  from 
20  to  40,  to  60,  to  40,  and  to  20  A  within  100  min. 

Fig.  13  shows  the  values  of  stack  voltage  in  terms  of  cur¬ 
rent  change  as  a  sine  function  I-  31+30  sin(/7r/30)  where  the 
stack  current  curve  was  plotted  as  a  reference.  During  the 
current  change  from  1  to  61  A,  the  voltage  output  slew  in 
the  range  27-38  V  and  the  minimum  voltage  output  value 
was  attained  when  the  current  (power  output  curve  had  the 


Fig.  10.  Transient  plot  of  stack  temperature  and  voltage  during  the  load 
set-up  from  30  to  50  A. 


Fig.  11.  Exit  temperature  change  with  time/load  (0a, in  =  0c, in  =  1.0). 


Fig.  12.  Stack  temperature  and  voltage  change  (0a,in  =  0c, in  =  1.0). 


Fig.  13.  Stack  voltage  as  a  function  of  current  which  changes  with  time  as 
I-  31  +  30  sin(t7i/30). 


194 


X.  Yu  et  al.  /  Journal  of  Power  Sources  147  (2005)  184-195 


Fig.  14.  Stack  efficiency  and  power  as  a  function  of  current  which  changes 
with  time  as  1=  31+30  sin(t7r/30). 

same  pace)  was  at  its  maximum  value  due  to  voltage  ohmic 
loss. 

For  the  steady  case,  the  average  efficiency  was  around 
45-65%,  depending  on  the  voltage,  energy  loss  to  the  sur¬ 
rounding  and  the  stream  sensible  heat.  For  the  unsteady  case, 
when  the  current  approached  zero,  the  efficiency  approached 
its  peak  value  with  the  maximum  attained  at  close  to  0A; 
the  efficiency  then  quickly  reduced  as  the  current  increased 
as  observed  in  Fig.  14.  When  the  current  and  power  output 
had  the  same  phase,  therefore,  when  the  current  reached  its 
peak  in  each  period,  the  power  output  also  reached  its  maxi¬ 
mum  value,  however,  stack  efficiency  had  the  opposite  trend 
(Figs.  13  and  14). 

5.  Conclusions 

The  main  objective  of  the  present  study  was  to  investigate 
the  Ballard  PEM  fuel  cell  performance  in  terms  of  electro¬ 
chemical  characteristics,  and  thermal  and  water  management. 
A  simplified  mathematical  model  and  tool,  which  could  be 
used  to  evaluate  PEM  fuel  cell  stack  electrochemical  perfor¬ 
mance  as  well  as  water  and  thermal  management,  has  been 
proposed  and  implemented  to  simulate  a  Ballard  PEM  fuel 
cell  stack. 

The  model  and  tool  can  be  applied  in  PEM  fuel  call  de¬ 
sign  processes  to  assist  the  designer  in  achieving  an  optimal 
design.  Knowing  a  set  of  gas  feeding  and  stack  physical  con¬ 
ditions,  the  design  engineer  can  obtain  the  information  re¬ 
garding  reaction  products,  stack  power,  temperature,  system 
efficiency,  and  dynamic  characteristic  of  the  PEM  fuel  cell 
stack,  by  use  of  the  numerical  trial-and-error  method  devel¬ 
oped  in  the  present  study. 

For  the  studied  Ballard  PEM  fuel  cell  stack,  the  more 
the  water  supplied  to  anode  from  its  inlet  (0a, in=  1-0  and 
1.5),  the  higher  the  voltage,  and  usually  the  lower  the  an¬ 
ode  exit  temperature.  The  studied  Ballard  PEM  fuel  cell 


stack  usually  takes  about  30-40  min  to  reach  its  steady  oper¬ 
ation. 
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